Sujet2-Etape1

1 Etape 1 : Première manipulation du RPG

1.1 Je choisis un point sur la carte : https://www.google.fr/maps

=> clic droit sur la carte + clic gauche sur les coordonnées (copiées dans le presse-papier)

J’indique un rayon en mètres (< 15000) pour sélectionner les parcelles autour du point

Code

# coord_gmaps<- rstudioapi::showPrompt(title = "Collez les coordonnées", message = "coordonnées Gmaps", default = "")
coord_gmaps<-"43.447894436406216, 1.2886163291688764"
lat<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,1])
lon<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,2])
# rayon<-as.numeric(rstudioapi::showPrompt(title = "Rayon", message = "Rayon (en m)", default = ""))
rayon<-10000

# création d'une table sf «point» avec les coordonnées saisies
# transformation des coordonnées en syst de proj 2154 (Lambert II - Français) 
point<-data.frame(lon,lat,rayon) %>% 
  st_as_sf(coords = c("lon","lat"),crs = "EPSG:4326") %>%
  mutate(coord_pt_gps=st_as_text(geometry)) %>% 
  st_transform("EPSG:2154") %>% 
  st_sf() %>% clean_names() %>% 
  rename(geom=geometry)

# st_crs(point)
# st_geometry_type(point)
# plot(point)

La table des parcelles agricoles 2021 se trouve sur un serveur PostgreSQL/PostGIS.

1.2 Je me connecte au serveur PostgreSQL avec le mot de passe disponible sur Onyxia/Mes services/PostgreSQL/Read me

Code

# le mot de passe est stocké dans un secret Vault

# postgresql_password <- rstudioapi::askForPassword(prompt = "Entrez le password PostgreSQL")
# postgresql_password <- "1tfawt3nj7fgzo3w7cma"
# Sys.setenv("PASS_POSTGRESQL"="1tfawt3nj7fgzo3w7cma")
  
# Connection à PostgreSQL
 cnx <- dbConnect(PostgreSQL(),
                user = "projet-funathon",
                password = Sys.getenv("PASS_POSTGRESQL"),
               host = "postgresql-758156",
               dbname = "defaultdb",
                port = 5432,
                options="-c search_path=rpg,public") # specify what schema to connect to

1.3 je lance une requête SQL pour sélectionner les parcelles situées autour de mon point dans le rayon choisi.

Code

# suppression de la table «point» si elle existe
dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.point CASCADE;")

<PostgreSQLResult>

Code

# écriture de la table point dans la base PostGis
write_sf(point, cnx, append = T)

# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.point ADD CONSTRAINT point_pkey PRIMARY KEY(coord_pt_gps);")

<PostgreSQLResult>

Code

# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")

<PostgreSQLResult>

Code

# dbExecute(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")

# 3 - Exécution de la requête de découpage du RPG autour du point sur PostGis  -------

dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.parc_prox CASCADE;")

<PostgreSQLResult>

Code

dbSendQuery(cnx,"CREATE TABLE rpg.parc_prox AS
    SELECT row_number() OVER () AS row_id, p.coord_pt_gps, p.rayon, r.*  
    FROM rpg.point p, rpg.parcelles r 
    WHERE ST_DWithin(p.geom,r.geom,p.rayon);")

<PostgreSQLResult>

Code

# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.parc_prox ADD CONSTRAINT parc_prox_pk PRIMARY KEY(id_parcel);")

<PostgreSQLResult>

Code

# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX parc_prox_geom_idx ON rpg.parc_prox USING gist(geom);")

<PostgreSQLResult>

Code

# 4 - Téléchargement des parcelles proches sous R-------------------------------

# parc_prox<-read_sf(cnx,parc_prox)

parc_prox<-st_read(cnx, query="select * from rpg.parc_prox;")

1.4 J’affiche les parcelles autour de mon point avec une carte interactive leaflet

Code

#plot(st_geometry(parc_prox))
#plot(st_geometry(point),add = T,col = "red")

# Transformation de la projection car leaflet ne connait que le WGS 84
carte_parc_prox<- parc_prox %>% st_transform(4326)

# Marqueur du point
pt_mark<- point %>% st_transform(4326)

# ajout du libellé des cultures
carte_parc_prox_lib<-carte_parc_prox %>% 
  left_join(lib_cult %>% select(-code_groupe_culture),by=c("code_cultu"="code_culture")) 
#création d'un label ad hoc à afficher en surbrillance au passage de la souris sur la carte.
labels<-sprintf("<strong>id_parcel : </strong>%s<br/>
                <strong>Groupe culture : </strong>%s<br/>
                <strong>Culture : </strong>%s<br/>
                <strong>Surface (ha) : </strong>%s<br/>
                <strong>Département : </strong>%s<br/>
                <strong>Commune : </strong>%s<br/>",
                parc_prox$id_parcel,carte_parc_prox_lib$libelle_groupe_culture,
                carte_parc_prox_lib$libelle_culture,parc_prox$surf_parc,
                parc_prox$insee_dep,parc_prox$nom_com) %>% 
  lapply(htmltools::HTML)

# labels

# création d'une palette de couleurs associée au groupe de culture
factpal <- colorFactor("Paired", parc_prox$code_group)

carte_parc_prox_html <- leaflet(carte_parc_prox_lib) %>% 
   # addProviderTiles("Esri.WorldImagery") %>%
   # addTiles() %>% 
  addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%
  addPolygons( #fillColor="white",
               fillColor=~factpal(code_group),
               weight = 2,
               opacity = 1,
               color = "#ffd966",
               dashArray = "3",
               fillOpacity = 0.5,
               highlight = highlightOptions(
                 weight = 5,
                 color = "#A40000",
                 dashArray = "",
                 fillOpacity = 0.0,
                 bringToFront = TRUE),
               label = labels,
               labelOptions = labelOptions(
                 style = list("font-weight" = "normal", padding = "3px 8px"),
                 textsize = "15px",
                 direction = "auto",
                 encoding="UTF-8")) %>% 
  addMarkers(data=pt_mark,~lon,~lat, popup = ~coord_pt_gps, label = ~coord_pt_gps)

Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Paired is 12
Returning the palette you asked for with that many colors

Code

# leaflet(pt_mark) %>% addTiles() %>% addMarkers(~lon,~lat)

carte_parc_prox_html

Code

# saveWidget(widget = carte_parc_prox_html, file = "carte_parc_prox_html.html")

1.5 Quelle est la structure des parcelles agricoles autour de mon point ?

Code

t1 <- parc_prox %>% st_drop_geometry() %>% count(code_group) %>% 
  add_tally(n) %>% 
  mutate(n_pct=round(100*n/nn,1)) %>% 
  select(-nn) %>% rename(n_parcelles=n) %>%
  # adorn_totals() %>% 
cbind(
  # comptage des surfaces
parc_prox %>% st_drop_geometry() %>% count(code_group,wt=surf_parc) %>% 
      add_tally(n) %>% 
      mutate(surf_pct=round(100*n/nn,1)) %>%
      select(-nn) %>%  
      rename(surf_parc_ha=n) %>% select(surf_parc_ha, surf_pct) # %>% 
       # adorn_totals()
) %>% left_join(lib_group_cult,by=c("code_group"="code_groupe_culture")) %>% 
  select(code_group,libelle_groupe_culture,everything()) %>% 
  # arrange(as.numeric(code_group)) %>% 
  arrange(desc(surf_parc_ha)) %>% 
  adorn_totals() %>% 
  mutate(taille_moy_parc=round(surf_parc_ha/n_parcelles,1))

Storing counts in `nn`, as `n` already present in input
Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.

Code

t1  %>% 
  setNames(c("code","groupe de cultures","nombre de parcelles","(%)","surface (ha)","surface (%)","taille moyenne (ha)")) %>% 
  kable(
    format="html",
    caption="<span style='font-size:medium'>Groupes de cultures <strong>locales</strong> par surfaces décroissantes</span>",
    format.args = list(decimal.mark = ",", big.mark = " "),
    booktabs = TRUE) %>%
  kable_styling(font_size = 15) %>% 
  gsub("font-size: initial !important;",
       "font-size: 20pt !important;",.)%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  row_spec(nrow(t1), bold = T, color = "white", background = "grey")

Structure des cultures au niveau local

code	groupe de cultures	nombre de parcelles	(%)	surface (ha)	surface (%)	taille moyenne (ha)
1	Blé tendre	631	10,9	3 044,52	17,6	4,8
2	Maïs grain et ensilage	368	6,4	2 922,39	16,9	7,9
4	Autres céréales	404	7,0	1 826,24	10,6	4,5
11	Gel (surfaces gelées sans production)	1 437	24,9	1 494,91	8,7	1,0
6	Tournesol	308	5,3	1 430,40	8,3	4,6
18	Prairies permanentes	573	9,9	1 260,24	7,3	2,2
7	Autres oléagineux	184	3,2	1 237,37	7,2	6,7
3	Orge	214	3,7	947,61	5,5	4,4
16	Fourrage	234	4,1	804,31	4,7	3,4
19	Prairies temporaires	302	5,2	802,74	4,6	2,7
8	Protéagineux	116	2,0	538,13	3,1	4,6
5	Colza	106	1,8	468,59	2,7	4,4
28	Divers	754	13,1	184,22	1,1	0,2
15	Légumineuses à grains	34	0,6	147,28	0,9	4,3
25	Légumes ou fleurs	38	0,7	46,95	0,3	1,2
21	Vignes	21	0,4	37,31	0,2	1,8
17	Estives et landes	16	0,3	33,19	0,2	2,1
24	Autres cultures industrielles	23	0,4	32,70	0,2	1,4
20	Vergers	7	0,1	4,76	0,0	0,7
Total	-	5 770	100,0	17 263,86	100,1	3,0

Code

# rm(t1)

1.6 Comparaison avec la répartition des cultures au niveau départemental et national

Code

# jointure spatiale avec la couche département pour récupérer le département où tombe le point

# ouvrir la couche des communes (à convertir en Lambert II epsg 2154)
dep <- s3read_using(
    FUN = sf::read_sf,
    layer = "departement",
    object = "2023/sujet2/diffusion/ign/adminexpress_cog_simpl_000_2023.gpkg",
    bucket = "projet-funathon",
    opts = list("region" = "")) %>% 
  st_transform(2154)

# jointure
df<-point %>% st_join(dep) %>% st_drop_geometry() %>% select(insee_dep)
dep_pt<-df[1,1]
rm(df)

# sinon sélection du département regroupant la plus grande surface agricole
# cas où le cercle recouvre plusieurs départements
# df<-parc_prox %>% st_drop_geometry() %>% 
#   count(insee_dep,wt=surf_parc) %>% 
#   arrange(desc(n))
# dep_pt<-df[1,1]
# rm(df)

# calcul des % surfaces autour du point
stat_pt <- parc_prox %>% st_drop_geometry() %>% 
  count(code_group,wt=surf_parc) %>% add_tally(n) %>% 
  mutate(pct_surf_local=round(100*n/nn,1)) %>%
  select(code_group, pct_surf_local) 

Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.

Code

# récup des % surfaces départementales
stat_dep_pt<-s3read_using(
  FUN=readr::read_rds,
  object = "2023/sujet2/diffusion/resultats/stat_group_cult_by_dep.rds",
  bucket = "projet-funathon",
  opts = list("region" = "")) %>% 
  filter(insee_dep %in% dep_pt) %>% 
  select(insee_dep,code_group,libelle_groupe_culture,pct_surf) %>% 
  rename(pct_surf_dep = pct_surf)

# récup des % surfaces nationales
stat_fm<-s3read_using(
  FUN=readr::read_csv,
  object = "2023/sujet2/diffusion/resultats/stat_group_cult_fm.csv",
  col_types = cols(code_group = col_character()),
  bucket = "projet-funathon",
  opts = list("region" = "")) %>% 
 select(code_group,libelle_groupe_culture,pct_surf) %>% 
  rename(pct_surf_fm = pct_surf)

# appariement des stas locales, départementales, nationales
stat_compar<-stat_fm %>% 
  left_join(stat_dep_pt %>% select(code_group, pct_surf_dep),by="code_group") %>% 
  left_join(stat_pt,by="code_group") %>% 
  select(libelle_groupe_culture,pct_surf_local,pct_surf_dep,pct_surf_fm) %>% 
  arrange(desc(pct_surf_local)) %>% adorn_totals() 

stat_compar %>% 
  setNames(c("Groupe de cultures","surf. locales (%)","surf. départ. (%)","surf. France m. (%)")) %>% kable(
    format="html",
    caption="<span style='font-size:medium'>Comparaison des surfaces locales, départemenales et nationales</span>",
    format.args = list(decimal.mark = ",", big.mark = " "),
    booktabs = TRUE) %>%
  kable_styling(font_size = 15) %>% 
  gsub("font-size: initial !important;",
       "font-size: 20pt !important;",.)%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  row_spec(nrow(stat_compar), bold = T, color = "white", background = "grey")

Structure des cultures

Groupe de cultures	surf. locales (%)	surf. départ. (%)	surf. France m. (%)
Blé tendre	17,6	14,3	17,7
Maïs grain et ensilage	16,9	8,6	10,0
Autres céréales	10,6	12,4	3,9
Gel (surfaces gelées sans production)	8,7	4,4	1,6
Tournesol	8,3	12,6	2,5
Prairies permanentes	7,3	17,9	27,7
Autres oléagineux	7,2	3,1	0,7
Orge	5,5	3,5	6,2
Fourrage	4,7	5,9	3,7
Prairies temporaires	4,6	4,2	5,3
Protéagineux	3,1	1,8	1,2
Colza	2,7	1,8	3,5
Divers	1,1	1,3	1,3
Légumineuses à grains	0,9	0,3	0,2
Légumes ou fleurs	0,3	0,3	1,5
Estives et landes	0,2	6,9	8,1
Vignes	0,2	0,4	2,2
Autres cultures industrielles	0,2	0,1	1,8
Vergers	0,0	0,1	0,4
Plantes à fibres	NA	0,0	0,4
Riz	NA	NA	0,0
Fruits à coque	NA	0,0	0,2
Oliviers	NA	0,0	0,0
Total	100,1	99,9	100,1

1.7 Graphique de comparaison des cultures au niveau local, départemental et national

Code

# je sélectionne les 10 groupes de cultures les plus répandus au niveau local 
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>% 
  rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)

# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name) 

# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)  
tab_piv[is.na(tab_piv)] <- 0


# levels(as.factor(tab_piv$secteur))
# je réordonne les secteurs dans le "bon" ordre, avec factor  
tab_piv$secteur<-factor(tab_piv$secteur,levels=c("france","departement","local"))
tab_piv<-tab_piv %>% arrange(desc(secteur), desc(value))

# je réordonne les cultures par surface décroissante au niveau local, avec factor
x <- tab_piv %>% filter(secteur == "local") %>% arrange(value) %>% select(libelle_groupe_culture)
y <- pull(x, libelle_groupe_culture)

tab_piv$libelle_groupe_culture <- factor(tab_piv$libelle_groupe_culture, levels = y)

# ggplot => problème de tri : affichage des prairies avant le maïs ?    
p<-ggplot(tab_piv, aes(x =libelle_groupe_culture,
                       y = value, 
                       fill=factor(secteur, levels=c("france","departement","local")))) + 
  geom_col(position = "dodge") +
  labs(title="Surfaces comparées des 10 principales cultures locales, en %", x="Culture", y = "%", fill = "Secteur") +
  theme_classic()

# je flippe le graphique pour avoir des barres horizontales  
p+coord_flip()

différences de structure locale, departementale, nationale

1.8 Je teste un graphique par secteur avec facet_wrap

Code

# je sélectionne les 10 groupes de cultures les plus répandus au niveau local 
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>% 
  rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)

# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name) 

# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)  
tab_piv[is.na(tab_piv)] <- 0

# ggplot => problème de tri : affichage des prairies avant le maïs ?    
ggplot(tab_piv, aes(x = reorder(libelle_groupe_culture,+value), 
                       y = value)) + 
  geom_col(fill = "lightblue", colour = "black",position = "dodge") +
  labs(title="Surface par culture", x="Culture", y = "%", fill = "Secteur") +
  geom_text(aes(label = value), hjust = -0.3, size = 8/.pt, colour = "black") +
  theme_classic() + coord_flip() + 
  facet_wrap(~secteur,nrow=3,scales='free')

différences de structure locale, departementale, nationale

1.9 Pour aller plus loin, comparer la taille moyenne des parcelles selon les secteurs